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cell,  which  acts  as  a  sink  and  source  reservoir,  is  maintained  at  reaction  equilibrium  conditions  via  the  RxMC  algorithm.  The 
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while  using  some  elements  of  the  osmotic  molecular  dynamics  method,  and  so  simulates  conditions  that  directly  relate  to  real, 
open  systems.  The  accuracy  and  stability  of  the  method  is  assessed  by  considering  the  ammonia  synthesis  reaction  N2 
+3H2<^2NH3.  It  is  shown  to  be  a  viable  method  for  predicting  the  effects  of  nonideal  environments  on  the  dynamic  properties 
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A  molecular  simulation  method  to  study  the  dynamics  of  chemically  reacting  mixtures  is  presented.  The 
method  uses  a  combination  of  stochastic  and  dynamic  simulation  steps,  allowing  for  the  simulation  of  both 
thermodynamic  and  transport  properties.  The  method  couples  a  molecular  dynamics  simulation  cell  (termed 
dynamic  cell)  to  a  reaction  mixture  simulation  cell  (termed  control  cell)  that  is  formulated  upon  the  reaction 
ensemble  Monte  Carlo  (RxMC)  method,  hence  the  term  reaction  ensemble  molecular  dynamics.  Thermody¬ 
namic  and  transport  properties  are  calculated  in  the  dynamic  cell  by  using  a  constant-temperature  molecular 
dynamics  simulation  method.  RxMC  forward  and  reverse  reaction  steps  are  performed  in  the  control  cell  only, 
while  molecular  dynamics  steps  are  performed  in  both  the  dynamic  cell  and  the  control  cell.  The  control  cell, 
which  acts  as  a  sink  and  source  reservoir,  is  maintained  at  reaction  equilibrium  conditions  via  the  RxMC 
algorithm.  The  reaction  ensemble  molecular  dynamics  method  is  analogous  to  the  grand  canonical  ensemble 
molecular  dynamics  technique,  while  using  some  elements  of  the  osmotic  molecular  dynamics  method,  and  so 
simulates  conditions  that  directly  relate  to  real,  open  systems.  The  accuracy  and  stability  of  the  method  is 
assessed  by  considering  the  ammonia  synthesis  reaction  N2+3H2^2NH3.  It  is  shown  to  be  a  viable  method 
for  predicting  the  effects  of  nonideal  environments  on  the  dynamic  properties  (particularly  diffusion)  as  well  as 
reaction  equilibria  for  chemically  reacting  mixtures. 

DOT  10.1103/PhysRevE.70.061103  PACS  number(s):  82.60.-s,  02.70.Ns,  05.10.Ln,  82.33. -z 


I.  INTRODUCTION 

Predicting  and  understanding  the  physical  effects  of  non¬ 
ideal  conditions  (strong  intermolecular  forces,  nanostruc- 
tured  media,  etc.)  on  chemical  reaction  equilibria  is  critical 
in  many  fields  of  science  including  mixture  separation  and 
purification  using  porous  solids,  catalysis,  plasma  physics, 
and  shock  physics.  The  reaction  ensemble  Monte  Carlo 
method  (RxMC)  [1-3]  is  a  powerful  simulation  tool  for 
studying  reaction  mixtures  and  is  uniquely  capable  of  pre¬ 
dicting  shifts  of  reaction  equilibria  caused  by  such  highly 
nonideal  environments.  The  RxMC  method  requires  as  input 
only  the  intermolecular  potentials  and  the  ideal-gas  partition 
functions  for  the  reaction  species  that  are  present.  Further¬ 
more,  the  method  does  not  require  a  reactive-type  potential 
that  mimics  bond  breakage  and  formation.  Applications  of 
the  RxMC  simulation  method  include  reactive  systems  con¬ 
fined  in  porous  materials  [4-10],  reactions  of  plasmas 
[11,12],  reactions  in  supercritical  fluid  solvents  [10],  reac¬ 
tions  under  shock  [13],  and  still  others  [14-17].  Deviations 
from  the  ideal-gas  phase  reaction  equilibria  caused  by  non¬ 
ideal  conditions  can  be  determined  for  quantities  such  as  the 
fluid  density,  pressure,  and  species  concentrations. 

However,  dynamic  properties  such  as  diffusion  coeffi¬ 
cients  cannot  be  determined  using  the  RxMC  method.  In  this 
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work,  we  introduce  a  simulation  method  that  can  determine 
the  dynamic  properties  of  reaction  mixtures  at  equilibrium. 
The  method  is  akin  to  the  grand  canonical  molecular  dynam¬ 
ics  method  [18-20]  and  the  osmotic  molecular  dynamics 
method  [21,22].  These  methods  use  a  control  cell  to  maintain 
the  desired  chemical  potential  while  a  dynamic  cell  is  used  to 
determine  the  dynamic  properties.  For  the  method  introduced 
here,  termed  the  reaction  ensemble  molecular  dynamics 
(RxMD)  method,  the  control  cell  is  used  to  maintain  the 
system  at  the  reaction  equilibrium  conditions  [1-3].  RxMC 
forward  and  reverse  reaction  steps  are  performed  in  the  con¬ 
trol  cell  only,  while  constant-temperature  molecular  dynam¬ 
ics  steps  are  performed  in  both  the  dynamic  cell  and  the 
control  cell.  The  dynamic  cell  is  in  direct  contact  with  the 
control  cell  so  that  particles  are  able  to  move  freely  between 
the  cells.  Since  the  entire  simulation  box  (consisting  of  the 
control  cell  and  the  dynamic  cell)  is  in  thermodynamic  equi¬ 
librium  due  to  the  molecular  dynamics  steps,  reaction  equi¬ 
librium  conditions  are  established  in  the  dynamic  cell  as  a 
consequence  of  the  physical  contact  between  the  control  and 
dynamic  cell.  In  this  scenario,  fluid  properties  in  the  dynamic 
cell  are  unaffected  by  the  stochastic  reaction  steps  occurring 
in  the  control  cell.  Therefore,  dynamic  quantities  of  reaction 
mixtures  such  as  the  velocity  autocorrelation  functions  and 
the  diffusion  coefficients  can  be  accurately  determined  in  the 
dynamic  cell.  These  dynamic  properties  are  more  precisely 
dynamic  equilibrium  properties  since  they  describe  correla¬ 
tions  at  different  times  along  an  equilibrium  trajectory  [23]. 
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FIG.  1.  Schematic  of  the  reaction  ensemble 
molecular  dynamics  method.  The  model  reaction 
2A<^B  is  occurring.  Molecular  dynamics  par¬ 
ticle  displacement  steps  (  )  occur  in  all 

cells,  while  reaction  ensemble  Monte  Carlo  reac¬ 
tion  steps  (  ◄••••#.  )  occur  only  in  the  control 
ceU.  The  dashed  lines  denote  the  portion  of  the 
dynamic  cell  in  which  the  dynamic  properties  are 
calculated. 


It  is  important  to  note  further  that,  analogous  to  the  RxMC 
method,  the  RxMD  method  predicts  the  physical  effects  on 
reaction  equilibria,  as  opposed  to  predicting  chemical  effects. 
The  RxMD  method  does  not  provide  reaction  rate  informa¬ 
tion;  separate  methods  must  be  used  for  this  purpose  (for  a 
review  of  these  methods  see  Santiso  and  Gubbins  [24].  De¬ 
spite  these  limitations,  the  RxMD  method  can  provide 
unique  insight  into  the  molecular-level  dynamic  behavior  of 
a  wide  variety  of  reacting  systems.  In  particular,  it  enables 
the  simultaneous  study  of  the  influence  of  nonideality  and 
nanostructure  on  both  diffusion  and  reaction  equilibria. 

II.  METHODOLOGY 

The  reaction  ensemble  molecular  dynamics  method  pro¬ 
vides  insight  into  the  dynamic  phenomena  of  reaction  mix¬ 
tures  by  combining  a  stochastic  simulation  method  (RxMC) 
with  a  deterministic  method  (constant- T  MD).  The  RxMC 
method  can  be  performed  at  constant-volume  or  at  constant- 
pressure  conditions.  The  constant- volume  version  requires 
two  types  of  Monte  Carlo  moves:  (1)  particle  displacements 
and  (2)  forward  and  reverse  reaction  steps.  Multiple  reac¬ 
tions  can  be  simulated  simultaneously  by  including  the  for¬ 
ward  and  reverse  reaction  steps  for  each  reaction  in  order  to 
maintain  stoichiometry.  The  constant-pressure  version  of  the 
RxMC  method  requires  the  additional  step  of  fluctuating  the 
simulation  cell  volume  to  achieve  the  desired  pressure.  Fur¬ 
ther  details  of  the  RxMC  method  can  be  found  elsewhere 
[1-3].  In  this  work,  we  demonstrate  the  RxMD  method  at 
constant- volume  conditions;  extension  to  constant-pressure 
conditions  is  straightforward. 

The  reaction  ensemble  molecular  dynamics  method  is  a 
direct  extension  of  the  reaction  ensemble  Monte  Carlo 
method.  In  the  RxMD  method,  Monte  Carlo  displacements 
steps  are  replaced  by  molecular  dynamics  time  steps,  while 
forward  and  reverse  reaction  steps  are  still  performed  in  a 
Monte  Carlo  fashion.  Reaction  steps  in  the  RxMC  method 
require  particle  insertion,  particle  deletion,  and/or  particle 
identity  exchange  in  the  simulation  cell.  Such  conditions  are 
typically  not  suitable  for  the  molecular  dynamics  technique 
since  the  deterministic  pathway  of  the  particle  trajectories 
will  be  disrupted  by  such  events.  To  avoid  these  adverse 
effects  to  the  particle  trajectories,  the  control  cell  and  dy¬ 
namic  cell  simulation  setup  described  above  is  implemented. 
In  practice,  this  entails  positioning  one-half  of  the  control 
cell  on  each  side  of  the  dynamic  cell.  Such  an  arrangement 
allows  for  the  use  of  periodic  boundary  conditions  in  all 


directions.  Moreover,  total  momentum  is  inherently  con¬ 
served  due  to  symmetry  of  the  total  simulation  box.  A  sche¬ 
matic  of  the  RxMD  method  is  given  in  Fig.  1,  where  the 
model  reaction  2A<^B  is  occurring.  In  Fig.  1,  molecular 
dynamics  time  steps  are  performed  in  both  the  control  and 
dynamic  cells  while  reaction  steps  are  performed  in  the  con¬ 
trol  cell  only.  Since  the  control  cell  is  in  direct  contact  with 
the  dynamic  cell,  the  particles  are  able  to  move  freely  be¬ 
tween  both  cells.  To  minimize  interface  effects,  the  proper¬ 
ties  of  the  fluid  in  the  dynamic  cell  are  determined  from  an 
interior  portion  of  the  cell,  which  is  away  from  the  control- 
cell-dynamic-cell  interface;  i.e.,  at  any  time  step  only  par¬ 
ticles  that  reside  within  the  interior  portion  are  included  in 
the  properties  calculation.  The  dashed  lines  in  Fig.  1  are 
intended  to  illustrate  this  procedure.  Last,  particle  velocities 
for  newly  inserted  particles  are  drawn  from  a  Maxwell- 
Boltzmann  distribution  at  the  appropriate  temperature 
[23,25]. 

III.  VERIFICATION:  APPLICATION  TO  AMMONIA 
SYNTHESIS 

We  demonstrate  the  RxMD  method  using  the  ammonia 
synthesis  reaction  N2+3H2<=^2NH3  at  constant  volume.  Due 
to  the  illustrative  nature  of  this  work,  the  reacting  species  are 
modeled  as  simple  spherical  particles  interacting  through  an 
exponential-6  potential,  where  electrostatics  contributions 
are  ignored  [26].  The  potential  was  truncated  and  shifted  at 
1.05  nm.  The  potential  parameters  are  given  in  Table  I, 
where  the  Lorentz-Berthelot  combining  rules  [27]  were  used 
for  unlike-pair  interactions.  Molecular  partition  functions 
were  determined  using  JANAF  thermochemical  data  tables 
[28].  All  simulations  were  initiated  from  an  initial  mixture  of 
600  H2  molecules  and  200  N2  molecules.  For  the  state  con¬ 
ditions  simulated  here,  these  relatively  large  system  sizes 
result  in  large  simulation  cell  volumes  that  help  to  minimize 
any  interface  effects  caused  by  the  dynamic-cell-control-cell 
boundary.  The  interior  two-thirds  of  the  dynamic  cell  were 


TABLE  I.  Exponential-6  potential  parameters  [26]. 


Species 

^core 

e/ks  (K) 

a 

NH3 

1.12 

3.72 

244.9 

12.0 

N2 

1.03 

4.17 

97.1 

13.0 

H2 

1.25 

3.49 

30.4 

11.2 
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TABLE  11.  Comparison  of  properties  calculated  by  the  reaction  ensemble  molecular  dynamics  method  (RxMD)  with  the  reaction 
ensemble  Monte  Carlo  (RxMC)  and  molecular  dynamics  (MD)  methods  at  various  temperatures.^ 


Method 

^conf 

(J/mol) 

P 

[MPa] 

density 

[g/cm^] 

x(N2) 

x(H2) 

(10^  m^/s) 

RxMC 

-3.92 

0.02976 

300.0  K 

0.0001916 

0.8823 

0.0294i 

0.0882i 

NA" 

RxMD  dynamic  cell 

-3.56 

0.03034 

0.0002028 

0.8854 

0.029 19 

0.08526 

552. 14 

control  cell 

-4.I4 

0.03118 

0.0001 869 

0.8837 

O.O3OI5 

0.08672 

NA 

MD-NVT 

-3.84 

0.03044 

0.000191 

0.882 

0.0294 

0.0882 

552.85 

RxMC 

-383.7i9 

14.834 

T=  600.0  K 

0.03827 

0.526i 

O.II84 

0.3552 

NA 

RxMD  dynamic  cell 

-384.2n 

14.775 

0.03998 

0.5273 

0.1158 

0.3574 

3.5846 

control  cell 

-381.926 

14.866 

0.03879 

0.526i 

O.II63 

0.3585 

NA 

MD-NVT 

-382.6i3 

14.835 

0.0382 

0.526 

0.118 

0.355 

3.5788 

RxMC 

-203.338 

55.98 

900.0  K 

O.O6687 

0.178i 

0.2052 

O.6I63 

NA 

RxMD  dynamic  cell 

-201.627 

55.45 

0.06548 

O.I8I5 

0.1974 

0.6228 

2.78I9 

control  cell 

-202.1 15 

56.56 

0.06606 

0.1733 

0.2097 

0.6192 

NA 

MD-NVT 

-202.7i8 

55.74 

0.0668 

0.178 

0.205 

0.616 

2.7887 

^Reported  uncertainties  determined  from  block  averages  where,  for  example  -383. implies  -383.7+1.9  [23]. 
^x{i):  mole  fraction  of  species  where  N  is  the  number  of  molecules. 

^NA=not  applicable  to  method. 


used  to  calculate  the  fluid  properties.  Cubic  simulation  cells 
were  used  and  periodic  boundary  conditions  were  imposed  in 
all  three  directions  for  the  total  system.  All  reported  pres¬ 
sures  were  calculated  using  the  virial  expression  [25]. 

A  standard  NVT  molecular  dynamics  method  was  em¬ 
ployed  with  the  equations  of  motion  solved  using  the  Verlet 
leapfrog  algorithm  [25]  and  implementing  the  damped  force 
method  of  Brown  and  Clark  to  maintain  constant  tempera¬ 
ture  [29].  A  time  step  of  Ar=3.0  fs  was  used  for  a  total 
simulation  time  of  2  ns  following  an  equilibration  period. 
The  ratio  of  attempted  reaction  steps  per  molecular  dynamics 
steps  was  40:1;  this  choice  was  based  on  a  series  of  short 
simulation  runs  that  were  made  during  the  development  of 
the  method,  but  may  vary  for  different  systems. 

The  ammonia  synthesis  reaction  was  simulated  at  a  series 
of  temperatures  under  known  gas-phase  conditions  [4].  Com¬ 
parisons  of  the  quantities  calculated  by  the  RxMD  method 
are  made  with  two  different  methods.  First,  the  thermody¬ 
namic  quantities  calculated  from  the  RxMD  approach  (e.g., 
density,  conflgurational  energy,  pressure,  and  species  concen¬ 
trations)  are  compared  to  quantities  calculated  by  the  RxMC 
approach.  Second,  the  dynamic  quantities  of  the  RxMD 
method  (e.g.,  velocity  autocorrelation  functions)  are  com¬ 
pared  to  quantities  determined  from  an  NVT  molecular  dy¬ 
namics  mixture  simulation. 

For  a  comparison  of  the  thermodynamic  properties,  a 
constant-volume  RxMC  simulation  was  performed  at  the 
same  temperature  and  initial  mixture  as  the  RxMD  simula¬ 
tion.  Table  II  presents  a  comparison  of  results  for  the  RxMD 
and  RxMC  approaches,  where  quantities  calculated  in  both 
the  dynamic  cell  and  the  control  cell  are  given.  As  expected, 
quantities  in  the  dynamic  cell  and  control  cell  are  the  same 


within  statistical  uncertainty,  ensuring  that  the  RxMD  has 
equilibrated  properly.  The  quantities  calculated  using  the 
RxMD  and  the  RxMC  approaches  also  agree  within  statisti¬ 
cal  uncertainty  for  all  the  temperatures  considered. 

For  a  comparison  of  the  dynamic  properties,  NVT  molecu¬ 
lar  dynamics  simulations  were  performed  at  the  temperature, 
volume,  and  species  concentrations  that  were  used  or  deter¬ 
mined  in  the  RxMC  simulations.  Self-diffusion  coefficients 
(D)  for  each  species  were  determined  from  the  velocity  au¬ 
tocorrelation  functions  [23]  in  both  the  A/^VY-MD  simulations 
and  in  the  dynamic  cell  of  the  RxMD  simulations.  Compari¬ 
sons  of  the  self-diffusion  coefficients  determined  by  both 
methods  for  NH3  are  given  in  Table  IT  Excellent  agreement 
is  evident  between  the  two  approaches,  with  similar  agree¬ 
ment  found  for  the  self-diffusion  coefficients  of  N2  and  H2. 

IV.  DISCUSSION 

We  have  presented  a  simulation  tool,  the  reaction  en¬ 
semble  molecular  dynamics  method,  to  study  the  dynamics 
of  chemically  reacting  systems.  The  RxMD  method  is  ca¬ 
pable  of  predicting  the  physical  effects  of  nonideal  environ¬ 
ments  on  dynamic  properties  of  the  equilibrium  mixture.  The 
RxMD  method  combines  the  reaction  ensemble  Monte  Carlo 
method  with  the  NVT  molecular  dynamics  technique,  result¬ 
ing  in  a  simulation  method  that  mimics  real,  open  systems 
for  many  experimental  situations.  The  method  as  presented 
here  effectively  combines  two  simulation  runs  into  a  single 
simulation,  since  the  same  result  could  be  achieved  by  run¬ 
ning  an  RxMC  simulation  followed  by  a  mixture  AT^Y-MD 
simulation  at  the  conditions  determined  from  the  RxMC  re¬ 
sult.  The  value  of  the  RxMD  method  lies  in  its  potential 
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application  to  other  scenarios — e.g.,  where  the  dynamic  cell 
is  modeled  as  a  porous  solid  and  two  control  cells  at  different 
thermodynamic  conditions  are  used.  Such  an  arrangement 
mimics  combined  reaction  and  adsorption  phenomena  rel¬ 
evant  to  membrane  reactors  and  fuel  cells,  as  well  as  to  vari¬ 
ous  possible  nanochemical  devices  [30]. 
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